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ABSTRACT 

We present the first study of the dynamical evolution of an isolated star cluster that combines 
a significant population of primordial binaries with the presence of a central black hole. We 
use equal-mass direct N-body simulations, with N ranging from 4096 to 16384 and a pri- 
mordial binary ratio of — 10%; the black hole mass is about one percent of the total mass 
of the cluster. The evolution of the binary population is strongly influenced by the presence 
of the black hole, which gives the cluster a large core with a central density cusp. Starting 
from a variety of initial conditions (Plummer and King models), we first encounter a phase, 
that last approximately 10 half-mass relaxation times, in which binaries are disrupted faster 
compared to analogous simulations without a black hole. Subsequently, however, binary dis- 
ruption slows down significantly, due to the large core size. The dynamical interplay between 
the primordial binaries and the black hole thus introduces new features with respect to the 
scenarios investigated so far, where the influence of the black hole and of the binaries have 
been considered separately. A large core to half mass radius ratio appears to be a promising 
indirect evidence for the presence of a intermediate-mass black hole in old globular clusters. 

Key words: stellar dynamics — globular clusters: general — methods: n-body simulations 
— binaries: general 



1 INTRODUCTION 

Over the last few years some tantalizing, but yet far from conclu- 
sive, evidence has been accumulating in support of the idea that 
some star clusters could harbor a central black hole (BH) with a 
mass of the order of 1O 3 M or more. Detection of such an inter- 
mediate mass black hole (IMBH) has been claimed for A/15 and 
Gl iGerssen et all2003tlGebhardt et all2002Ll2005l) . However, al- 
ternative dynamical mod els without a central BH have been pro- 
posed for these clusters iBaumgardt et al. 2003a b). Interestingly, 
the visual appearance of globulars containing an IMBH is not that 
of a so-called core-collapsed clust er, but ra ther that of a cluster with 
a still sizable core ( Baumga rdt et all2005t) . 

IMBHs present a high theoretical and observational inter- 
est as these could be potential ultra-luminous X-ray sources and 
even emit gravitational waves, detectable by the next generation 
of gravitational wave detectors, as a result of close interactions 
with stars. However, despite this interest and the fact that theoret- 
ical studies of BHs in stellar systems started more than 30 years 
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ago (e.g., see |Peebles|[l972l iBahcall & Wol3ll976l) . detailed di- 
rect N-body simulations to study the dynamics of an idealized 
model with sing le stars and a central BH have been performed 
only recently iBaumgardt et alj|2004albt see also the studies on 
the formation of IMBHs by runaway mergers of massive stars by 
IPortegies Zwart et all2004 . 

One ingredient that complicates N-body simulations of glob- 
ular clusters is the presence of primordial binaries. Often, these 
are neglected in large simulations despite the increasing observa- 
tional e vidence that many stars in a globular cluster have a com- 
panion <Hut et alll992tfAlbrow et all200lUBellazzini et all2002t 
IZhao & Bai!vnTl200^r This frequent neglect is due to the dra- 
matic increase in computational resources required in a simula- 
tion where the local dynamical timescale may be many orders of 
magnitude smaller than the global relaxation timescale (hard bi- 
naries have an orbital period of a few hours, while the half-mass 
relaxation time can be up to a few billion years). The study of the 
dynamics in the presence of primordial binaries has been mainly 
limited to Fokker-Planck or Monte Carlo approaches iGao et all 
Il99lt iGiersz & Spurzemll200Ct iFregeau et 311 12003 ) and to direct 
simulations with rathe r modest particle numbers, from N ~ 10 3 
( McM illan et ailll990l: iMcMillan & Hu3ll994 iHeggie & Aarsethl 
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1992) t o recent higher resolution simulations , with TV up to 
16384 jHeggie. Trenti & Hutl 120061 : iTrenti. Heggie & Hutl l2006l 
hereafter Paperl and Paperll respectively). The results of these 
direct simulations have to be extrapolated to a realistic number 
of particles (e.g., by means of the IVesperini & Chemofj 1 19941) 
model) before being applied to the interpretation of the evolu- 
tion of typical globular clusters, that have 10 5 — 10 6 stars. Some 
realistic simulations, including primordial binaries are available 
IPortegies Zwart & Mc Millan 2004), but these are limited only to 
the first stage of the life of young dense clusters. In the case of 
open clusters, M67 has been modeled in a 3 6,000-body simulation 
running for several Gigayears iHurlev e t al. 2005). 

The presence of either an IMBH or a significant population 
of primordial binaries leads to an early release of abundant energy, 
inhibiting the development of a deep core collapse and hence of 
the onset of gravothermal oscillations. An IMBH can generate en- 
ergy by swallowing or tightly binding stars deep in its potential 
well [note also that energy can be generated through encounters 
of stars in the density cusp that is formed around the BH], while 
primordial binaries can generate energy by rapidly increasing bi- 
nary binding energy through three and four body encounters. En- 
ergy thus generated in the core of the system fuels the expansion 
of the half-mass radius, leading to a self-similar expansion of the 
entire system iHenon 1965). 

In this work we present the first direct simulations (with a 
number of equal-mass particles up to TV = 16384) of the evolu- 
tion of star clusters with both a significant fraction of primordial 
binaries (10 %) as well as an IMBH with a mass of 0.8 — 3 % of 
the total mass of the system. We address the following questions. 
Under the combined effect of the BH and of the binaries, what is the 
equilibrium size of the cluster core? How is the binary population 
affected by the presence of the central IMBH? Which physical pro- 
cesses dominate in the core? In principle two competing effects are 
possible: the BH may either enhance the disruption rate of binaries 
both indirectly, du e to the crea tion of a density cusp, and directly, 
by tidal stripping iPfahl 2005), or it may reduce the probability 
of interactions between binaries and sin gles, by producing a low 
stellar density in a relatively large core iBaumgardt et alj|2004al) . 
Which process is dominant is of fundamental importance. These 
questions are addressed in the next sections. 



2 NUMERICAL SIMULATIONS: SETUP 

Like in Paperl and Paperll, the simulati ons (see Tab le have 
been performed using the NBODY-6 code l Aarseth 2003). For this 
project NBODY6 has been modified with the kind help of Dr. 
Aarseth to ensure a more efficient and accurate treatment of the 
dynamics around the BH by fine tuning the parameters controlling 
chain regularization. 



2.1 Units 

All our results are presented 
(Heggie & Mathieu 1986) in which 



using standard units 



G = M 



-4E T 



where G is the gravitational constant, M is the total mass, and Et 
is the total energy of the system of bound objects. In other words, 
Et does not include the internal binding energy of the binaries, 
only the kinetic energy of their center-of-mass motion and the po- 
tential energy contribution where each binary is considered to be a 



point mass. We denote the corresponding unit of time by 

t d = GM 5/2 /(-4fi T ) 3/2 = 1. 

For the relaxation time, we u se the following expression I Spitzei 
ll987tlGiersz & Heggie! 1991) : 



trh = 



0.138TW 



3/2 



IGM In (0.11 TV) 

where r;, is the half-mass radius and TV denotes the number of orig- 
inal objects (binaries + singles). 

The core radius of the system is defined as the density aver- 
aged radius: 



(1) 



where r n is the distance of the n-th star from the density center, and 
the density p n around each particle is computed from the distance 
to the fifth nearest neighbor (Casertano & Hut 1985). A discussion 
o f other pos sible definitions for the core radius is extensively given 
in lCasertano & HutNl985l) : see also Paperll. 

Like in Paperl and Paperll, the internal energy of the binaries 
is measured in units of kT, where (3/2)fcT is the initial mean ki- 
netic energy per particle of the system (the binaries being replaced 
by their barycenter). 

As we are considering simulations with purely classical grav- 
itational interactions only, all our results are intrinsically scale-free 
and do not depend on the specific choices in physical units for the 
total mass and for the scale radius of the system modeled. 

For reference in physical units, a globular cluster (described 
by a Plummer model) with TV = 3 • 10 5 stars, a total mass of 
M = 3 • 10 5 Mq and a half-mass radius of 4 pc has a half-mass 
relaxation time t r h ~ 8.5 ■ 10 8 yr; in this cluster a binary, formed 
by equal mass stars each of mass 1 Mq, with binding energy of 
1 kT has a semi-major axis of m 10 All . 



2.2 Initial conditions 

The models considered in this paper are isolated with stars of equal 
mass m, plus a BH, introduced as a massive star, with mass (tush) 
in the range 0.8 — 3 % of the total mass of the system. The initial 
distribution is either a Plummer model (entries "PL" in Tab.0 or a 
King model (entries "K" in Tab.0 with concentration index Wo = 
3, 5, 7, 11. Our standard models have a primordial binary ratio of 
0-10 %. 

We define the primordial binary fraction as: 



/ = n b /(n B + rib) 



(2) 



with n s and rib being the number of singles and binaries respec- 
tively. 

We have considered runs with TV in the range 4096 — 16384. 
Note that TV denotes the number of original objects, i.e. TV — n s + 
rib, while the total number of stars is Ntot = n s + 2rib. when we 
discuss a run with TV = 8192 and 10% primordial binaries we are 
dealing with 9011 stars. 

The initial binding energy distribution for the binaries is in 
the range 5 to 680 kT, flat in logarithmic scale, similar to the 
initial conditions considered in the studies focused on the evolu- 
tion of star clusters wit h primordial binaries, but without an IMBH 
iGao et alll99lHFregeau et all20 03: Paperl; Paperll) 1 . This form 

1 Note that with our definition of kT we are considering the same range 
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Table 1. Summary of N-body runs 



ID ICtype N m BH /M f 

Pa Plummer 8192 0.014 0.1 

Pb Plummer 8192 0.014 

Pc Plummer 8192 0.025 0.1 

Pd Plummer 8192 0.1 

Pe Plummer 8192 

Pf Plummer 16384 0.1 

Pg Plummer 16384 0.008 0.1 

K5al King Wo = 5 8192 0.014 0.1 

K5a2 King Wo = 5 8192 0.014 0.1 

K5bl King Wo = 5 8192 0.014 

K5b2 King W = 5 8192 0.014 

K7al King Wo = 7 8192 0.014 0.1 

K7a2 King W = 7 8192 0.014 0.1 

K7bl King Wo = 7 8192 0.014 

K7b2 King W = 7 8192 0.014 

K3a King Wo = 3 4096 0.03 

K7a King W = 7 4096 0.03 

Klla King Wo = 11 4096 0.03 

Initial conditions for our set of runs. Columns (left to right): identification 
mark for the simulation, type of initial density profile, number of particles, 
black-hole mass, fraction of primordial binaries. 

for the binding energy distribution is suggested by the observed 
properties of binaries in star clusters iHut et all 19921) and the upper 
and lower limits are physically motivated: very soft binaries, i.e. bi- 
naries with binding energies below a few kT, would be quickly de- 
stroyed (i.e. on a timescale below one relaxation time) due to three 
and four body encounters in the core of the the system; binaries 
harder than a few hundreds kT are, on the other hand, dynamically 
inert. 

To initialize the simulation in a situation of approximate dy- 
namical equilibrium in the presence of the BH we first generate 
our chosen initial configuration (Plummer or King model) made 
of Ntot single stars only and we add the IMBH at rest at the cen- 
ter of the system. We then scale the velocities of the particles to 
reach virial equilibrium and let the system evolve for 5 half-mass 
crossing times td- In this first phase the system settles down in dy- 
namical equilibrium by experiencing mild virial oscillations (with 
amplitude below 6% of the equilibrium virial ratio) that are damped 
out in 2 — 3 td- At this stage we select randomly rib single stars that 
are eliminated from the simulation. Other rib single stars are ran- 
domly chosen to become binaries: for these a companion star is 
added with a semi-major axis (a) extracted from the assumed flat 
distribution in binary binding energy (that is translated in a dis- 
tribution in a proportional to 1 /a). The initial eccentricity for the 
binaries is chosen from a thermal distribution. At this stage the ini- 
tialization is complete with the system in approximate dynamical 
equilibrium. 

Our initial configurations have a limited choice of initial 
density profiles (Plummer model and King profiles with Wo = 
3, 5, 7, 11). This does however not influence or bias significantly 
the long-term, collisional evolution of the system. In fact, as we 
have shown in PaperII(see also Sec. 13.21 . the memory of initial 
conditions is erased within a few half mass relaxation times. Runs 

of lGao et alj ll99ll) . that reports the binary binding energy in units of kTc 
([3fcT c ; 400fcT c ]), with a temperature T c defined in terms of the central 
velocity dispersion (see Sec. 3.1 in Paperl) 



starting from concentrated King models initially presents a core 
expansion, while runs starting with shallower cores have an initial 
contraction, with both processes leading to a common dynamical 
configuration that can be interpreted in terms of balance between 
energy production in the core and energy dissipation due to half 
mass radius expansion (see alsofVesoerini & Chernofil 19941) . This 
argument also applies to justify our choice in the initialization pro- 
cedure to distribute, like in Paperll, primordial binaries so as to 
exactly trace the distribution of single particles, without exploring 
the possibility of a initial mass segregation. 

The evolution of the system is then followed up to t ~ 
25 t r h(0) that corresponds to about 20 Gyga year in physical units 
for a globular cluster with N = 3 ■ 10 s , M = 3 • 10 5 Mq and 
rh — 4 pc at t = 0. For t r h(0) we mean t r h computed at t — 0. 
For our initial conditions with N = 8192, t r h(0) ~ 112^. 

At N — 8192 the CPU time required for this typical run is 
of the order of 10 days on a state of the art pc linux workstation. 
As here we are interested in the collisional evolution of the sys- 
tem our results are expressed in terms of the relaxation time. This 
means that the computational complexity for direct N-body simu- 
lations, if lasting for a constant number of relaxation times, scales 
as N 3 /ln(0.11N): the factor iV 2 comes from the direct gravita- 
tional interaction, while the additional N/log (0.11 N) is due to the 
increase of the two body relaxation time with respect to the dynam- 
ical time jHeeaie & Hul2003l Chapter 3). 

2.3 BH properties and N scaling 

The mass of the black hole in our runs is in the range 
[0.8%; 3%] of the total mass of the system, a value that is 
slightly higher tha n the expected mass for a n IMBH following the 
mBH — cr relation jFerrarese & M erritt 2000 ; Gebhard t et alj200ol ; 
iTremaine et al|200 2j) app lied to a typical globular cluster. E.g. us- 
ing |Tremame_etal. (2002) with a velocity dispersion of 10km/ s 
we obtain a tjibh ~ 10 s Mq, which is about 0.3% to 1% of the 
total mass of a typical cluster. Our choice for the ratio ttibh/M 
is motivated by the limits that the current combination of hardware 
and software set on the number of particles that can be employed 
in direct simulations of gravitational systems with the combined 
presence of a black hole and primordial binaries. 

The number of particles that we can include in one simulation 
is between one and two orders of magnitude less than the typical 
number of stars in a globular cluster, so we are forced to introduce a 
scaling for the mass of the black hole. In principle two possibilities 
are available: we can keep constant either the ratio of the black hole 
mass to the stellar mass or the black hole mass to the total mass of 
the cluster. The first choice has the advantage of keeping unchanged 
the local interaction of stars with the black hole, while the second 
allows for a better comparison with observed star clusters. With 
respect to the properties of the interaction between the black hole 
and the primordial binary population it is not clear whether there 
is a preferred scaling. In fact, a black hole mass which is too big 
with respect to the total mass of the cluster may bias the results 
of the simulation, since the black hole may well be more massive 
than the core of the system, which usually contains a few percent 
of the total mass (see[Baumgardt et al. 2004aj). But similarly it can 
be argued that also a low ratio of the black hole to stellar mass in- 
fluences the results of the simulation, as it decreases the efficiency 
of tidal disruption of binaries around the black hole cusp, which 
is proportional to rriBH/rn- In addition a black hole with a ra- 
tio rriBH/m « 50, i.e. a 1 % mass black hole in a system with 
5 10 3 particles, has a small, but non negligible, probability of be- 
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ing ejected from the core as a result of a strong interaction with a 
binary. E.g. this happened during some test runs with tyibh = 0.01 
and TV = 4096 on a timescale of about t = 20 t r h(0), so that we 
were forced to stop these simulations at the moment of the ejection. 

Finally here we consider the black hole as a massive particle, 
and no mass accretion due to tidal disruption of stars is taken into 
account. However, this is not expected to bias the results signif- 
icantly: in fact (1) the mass accretion, whose rate is proportional 
to the ratio (vtl bh /m) 61 ^ 27 , is not ver y important in runs with 
TV < 10 4 (see iBaumgardt et aljE004l Table 1); in addition (2) 
stars passing close to the black hole may be captured, so that they 
become inert with respect to the other stars (effectively behaving 
like as if they had been accreted onto the black hole) or ejected from 
the system at extremely high speed (above 100 times the mean ve- 
locity dispersion), so that also in this case the star is removed from 
the system (with the difference that its mass is not added to that of 
the BH). 



3 GLOBAL EVOLUTION 

After the initialization described in Sec. 12.21 our simulations start 
from a condition of dynamical equilibrium: the evolution of the 
system is due to collisional effects and happens on a t r h timescale. 
This means in particular that the system is in virial equilibrium 
(2K/U — 1, where K is the kinetic energy and U is the poten- 
tial energy of the system). Only occasionally there are fluctuations 
of 2K/U at the order of one percent, when a strong interaction 
between a hierarchical system (e.g. a binary interacting with the 
IMBH) takes place leading to the consequent ejection of a high 
velocity particle from the system. The virial ratio returns to the 
equilibrium value on a dynamical time-scale. 

The large scale structure evolution of the cluster is dominated 
by the heating related to the presence of the BH. Starting from a 
Plummer model only the inner regions of the system experience a 
mild collapse on a time-scale that, depending on the mass of the 
BH, is of the order of a few t r j, (depicted in Fig.Q. Inside the core 
radius r c a cusp in the density profile is formed within the sphere 
of infl uence rj (with r, ~ 15 r c mBH /M, see Baumga rdt et all 
2004a) of the BH, with a profile proportional t o ~ 1/r 1 7 (see 
fig-GJ and thus similar to the 1 /r 1,75 measured by Baumga rdt et alJ 
l2004ap. For tubh = 0.014M the influence radius is approxi- 
mately 0.2r c . By definition, the stellar mass within this radius is 
comparable to that of the BH, and thus around one percent of the 
total mass of the cluster. 

The limited number of particles that we employ does not al- 
low us to significantly characterize the anisotropy profile of the sys- 
tem, that appears to be consistent with quasi-isotropy, except for the 
outer parts of the system, where there is a mild excess of radial or- 
bits. As the presence of a population of binaries does not influence, 
to a first approximation, the evolution of the anisotropy profile, we 
refer to the discussion in lBaumgardt et ail l2004at) (see in particu- 
lar their Fig. 8), where the authors take advantage of the absence of 
primordial binaries and run simulations with O(10 5 ) particles. 

3.1 Plummer Models 

For a Plummer model with TV = 8192 and / = 10%, in case of 
m,BH = 0.014M, the core radius r c is reduced in w 4 t r h(0) 
from the initial value of 0.4 to 0.3 (in natural units, see Sec. l2.1> . 
This is to be compared with a value of m 0.1 reached without the 
BH (Paperl; see also Tab.|2j; when a BH (with the same mass) but 
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Figure 1. Lagrangian radii (enclosing from bottom to top 2,5,10,20,50,80 
and 90 % of the total mass of the stars of the system) for a run with N=8 192 
starting from a Plummer profile with / = 10% and rriBH = 0.014. 



1 r 




Figure 2. Density profile (in natural units) vs radius (in units of the half 
mass radius r% = 1.44) at time t = 13 t r h(0) for a run with N=16384 
starting from a Plummer profile with / = 10% and tubh = 0.014. The 
plot shows the stellar mass only. The BH contribution to the plot would 
be that of a delta function located at r = 0. The large core size, basically 
that of a Plummer model, is an indirect consequence of the presence of the 
central IMBH. 
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Table 2. Global properties of Plummer runs with JV = 8192 



tubh/M 


/ 


(r c )cc/r c o 


r c /r h 
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0.29 
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0.1 


0.25 


0.09 


2.25 








0.02 
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In the first column we report the BH mass, in the second the fraction / of 
primordial binaries, in the third the core radius at the end of the initial core 

contraction phase (r c ) cc in units of the initial core radius r c rj, the fourth 
entry is the core to half mass radius ratio during the self-similar expansion 
of the system, while the last entry (r^f ) is the value of the half mass radius 
at t = 24 t r y l (0) in units of the initial half mass radius. 



no primordial binaries are present r c goes down to « 0.28 . From 
these values it is apparent that the IMBH dominates the global evo- 
lution of the system, as there is very little difference in simulations 
with and without primordial binaries if an IMBH is present. 

After the first mild contraction, all Lagrangian radii (shown 
in Fig.0 start to expand steadily and a self-similar regime sets in, 
with the half-mass radius growing in proportion to ~ t 2 ' 3 , in agree- 
ment with the theoretical argument given by Henon 1 1965) and as 
found for runs with single stars only by Baumsardt et al. 1 2004a). 
The rate of expansion of the half-mass radi us is set by the amo unt 
of energy generated in the core IVesperini & Chernofj[l994l) . In 
our runs with binaries plus IMBH this is marginally bigger (by 
~ 10 % at 24 trfc(O)) than in runs with binaries only; without a 
BH the half-mass expansion starts only after core collapse, which 
takes ss 10 t r h(0) in that case (Paperl). Conversely, if a BH but no 
binaries are present, the expansion rate of r% is reduced by « 20 %. 
A summary of the properties of the Plummer runs is reported in Ta- 
ble|2| 

This picture does not appear to depend much on the number 
of particles used, except quantitatively for the known logarithmic 
dependence on A of the ratio r c /rh IVesperini & Chernoffll994T) . 
To validate the results of our series of A = 8192 runs we have 
performed a run with A = 16384 and an IMBH mass tubh = 
0.008 (thus keeping quasi-constant the ratio m/rriBH and going 
toward a BH mass closer to the one expected by the rriBH — & 
relation). The evolution of the system, reported in Fig. [5] is indeed 
very similar, both in terms of binary destruction rate and in terms 
of the Lagrangian radii evolution. 



3.2 King Models 

At fixed number of particles, IMBH mass and number of primor- 
dial binaries, there is no expectation that the long term evolution 
of the system depends much on the details of the initial conditions, 
as these are erased on a relaxation time-scale (see also Paperll). 
To verify this we have also started some of our runs (see Table0 
from an initial King profile, with different values of the concentra- 
tion parameter Wo- The evolution of the core radius with respect 
to the half mass radius is reported in Fig. [4] Both the Wo = 5 and 
the Wo = 7 models evolve toward r c /r^ w 0.3 on a time-scale 
of « 5 trh(0) (this is the same value reached in the simulations 
starting from a Plummer profile). Interestingly the Wo = 7 model 
presents an initial expansion of the core, as the initial central con- 
centration is so high that the energy released by three and four body 
encounters cannot be completely dissipated by expansion of the 
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Figure 3. Comparison between two simulations starting from a Plummer 
profile with / = 10%: one has N = 16384 with m BH = 0.008 (solid 
line), the second has N = 8192 with ra BH = 0.014 (dashed line). In 
the upper panel we depict the relative number of binaries and singles, in 
the lower panel selected lagrangian radii (enclosing 5, 10, 20, 50 and 80 % 
of the total mass in stars). The Lagrangian radii have been smoothed with 
a triangular window filter of width 1.5 t r f l {0). Doubling the number of 
particles does not introduce significant changes in our simulations. There 
is only a small contraction of the Lagrangian radii, when the BH mass is 
reduced to 0.008. 
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Figure 4. Core to half mass radius ratio for a series of simulations with 
N = 8192, tubh = 0.014 and 10% primordial binaries. Different ini- 
tial conditions (King profiles with Wo = 5,7) converge toward a common 
value for r c /r h . The results are presented by applying a triangular smooth- 
ing window of width 1.5t r h(0). 
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half mass radius. This behavior is similar to that observed in runs 
with primordial binaries only but with yet higher central concentra- 
tions iFregeau et al 120 03: Paperll). This picture is confirmed also 
in a different series of runs - starting from Wo = 3, 5, 7, 11 with 
itibh = 0.03 and no binaries -that exhibit a very similar behavior. 



4 PROPERTIES OF THE BINARY POPULATION 

The evolution of the number of binaries in our simulations is driven 
by two competing processes: the formation of new binaries and 
the destruction of existing ones. As we discussed in detail in Pa- 
perl (see in particular Sec. 2 in that paper), the probability of form- 
ing new binaries is greatly suppressed with respect to the proba- 
bility of destructing an existing binary. In fact, the leading forma- 
tion channel for a new binary from single stars is a three body en- 
counter, which is proportional to p 3 . p s is the number density of 
single stars and this density is expressed in physical units, appro- 
priate for encounters, such that p = 1 corresponds to inter-particle 
distances equal to the ninety-degree turnaround impact parameter 
1989). Therefore in this units we have p s <C 1. The leading 
destruction channels are instead given by binary-single and binary- 
binary encounters and therefore proportional to pbPs and pi, where 
Pb is the number density of binaries, in the same physical units con- 
sidered for p a . As during our simulations the core number density 
of binaries is similar to that of singles, these two last channel domi- 
nates over the three single stars encounters. This expectation is con- 
firmed empirically in our runs, where the number of dynamically 
formed binaries that can be found after a few relaxation times in 
the simulation (typically about 10 -3 of the number of single stars) 
is negligible not only with respect to the number of primordial bi- 
naries but also with respect to the number of disrupted binaries. 
Essentially, the dominant dynamical mechanism that changes the 
net number of binaries is that of binary disruption. Of course, the 
exchange of binary members during a binary-single encounter is 
also frequent, but this does not alter the net number of binaries in 
the cluster. 

As the presence of the BH dominates the global dynamics of 
the cluster, the evolution of the binaries presents some remarkable 
differences from the scenario where no central BH is present. We 
can distinguish two phases. 

The presence of a BH significantly accelerates the rate at 
which binaries are disrupted in the first few half-mass relaxation 
times (see Fig.[3J- The reason is that binaries that pass near the 
BH can be quickly destroyed. The disruption rate is approximately 
constant during the first 5 t r h(0), while without a BH disruption 
takes a while to get underway. The initial disruption rate depends 
on the mass of the BH: a more massive BH starts burning binaries 
at a higher rate, as can be seen from Fig. [5] 

For our Plummer run with N = 8192, tubh = 0.014M and 
/ = 10% we observe a binary depletion rate of « 0.15 binaries per 
initial half-mass crossing time (td), during the first 5 t r h (0) (here 

trh(0) = 116.1t d ). 

For this time interval we present in Fig. [6j the distribution of 
the radial positions at which binaries are destroyed. This is ob- 
tained by looking in NBODY6 at the position at which binaries are 
flagged by the code as too wide to continue meriting special nu- 
merical treat ment in terms of KS I Kustaanheimo & Stiefel 1965) 
solutions (see lAarsethll2003l for a detailed description of the algo- 
rithm). In short, a KS solution is terminated in case of an external 
perturbation at the level of 5% or higher of the internal gravitational 
force between the members of the binary (if the binary is not de- 
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Figure 5. Time-dependence of the number of single (upper set of lines) and 
binary stars (lower set of lines) expressed as a fraction of the initial values. 
The data refer to simulations with TV = 8192 and 10% primordial binaries. 
There are three choices for BH mass: no BH (dotted line), a central BH with 
rriBH = 0.014M (solid line), and m BH = 0.025M (dashed line). The 
unit of time is the initial half-mass relaxation time. 



stroyed by the perturbation, its special numerical treatment is even- 
tually restarted when the external force returns below the critical 
5% threshold). In the histogram of Fig.|6|there are 226 events reg- 
istered, but these include also close encounters with perturbers that 
do not change the net numbers of binaries (e.g. exchanges). The 
net number of destroyed binaries for the time interval considered 
is 80 (see Fig.|5J. Around half of the binaries are disrupted within 
the black hole sphere of influence n (100 out of 226 KS solutions 
terminated within r;). Interestingly only two of them happen to be 
so close to the BH tidal radius rt that the disruption may be consid- 
ered related to the tidal stripping force exerted by the BH (see left 
panel of Fig.|5J. Here the tidal radius is defined as the radius where 
the gravitational force between the members of the binary equals 
the tidal force exerted by the BH (see e.g., Spitz eill987l) . i.e. as 

r< = (m B H/m) 1/3 a, 

where a is the semi-major axis of the binary. With t his respect, the 
theoretical model recently proposed by Pfahl (2005) applied to our 
N = 8192, m B H = 0.014, / = 10% simulation would give 
(from his Eqs. 1 1-12) a tidal disruption rate of « 2 ■ 10 -2 binaries 
per initial half-mass crossing time, so that we observe less tidally 
stripped binaries than expected (w 10). 

The disruption of binaries by the BH is thus mainly due to 
an indirect effect: binaries venturing close to the BH, where the 
density is higher, are more likely to interact with a single star or 
with another binary before having the chance to being tidally de- 
stroyed by the IMBH. In fact in our simulations we often observe 
the presence of hierarchical systems within rj, with the BH playing 
the role of a perturber. Fig.|6|nicely summarizes this picture: bina- 
ries are disrupted at a significantly enhanced rate within n (right 
panel), but the minimum distance from the BH is the order of 10 rt 
(left panel). The bi-modal distribution in the left panel of Fig. |S| 
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also reflects the enhanced probability of close encounters and bi- 
nary disruption within the influence sphere of the BH, as can be 
seen from Fig.0 

After the initial transient phase a self-similar expansion sets 
in, where the average core density is much less (approximately by 
a factor 10) than it would be without the presence of a central com- 
pact object. In this second phase we observe (see Fig.[5J a reduced 
rate of binary disrupti on (as this rate is prop ortional to the square 
of the density, e.g. see IVesperini & Chemofi 19941) . so that the dif- 
ference between a simulation with and without a BH becomes re- 
markable. The turning point is around 10 t r h(0), when the number 
of surviving binaries for a simulation with a central object becomes 
greater than in the absence of a BH. Our Fig.|5]has been given in 
units of the initial half-mass relaxation time, but the picture remains 
qualitatively the same even if we consider a co-moving time coor- 
dinate to take into account the differences between the half-mass 
radii of a simulation with and without a BH. 

Interestingly the spatial distribution of binaries is also differ- 
ent from the case without a BH (see Figs.|8j. After w 20 t r h{0) the 
number of binaries within the 0.05 Lagrangian radius is much less 
for the simulation with a BH. This is probably due to the disruption 
of binaries that approach close to the BH. In the presence of a BH, 
binaries tend to be more concentrated between the 0.05 Lagrangian 
radius and the half-mass radius, while in the absence of a BH, bina- 
rie s can sink deeper into the central region of the cluster. As noted 
bv lBaumgardt et all j2004a) for simulations without primordial bi- 
naries, the mass segregation efficiency is different in presence of a 
central BH. 

If we compare the evolution of the binding energy distribu- 
tion of the binaries in a run with and without a central BH, we 
can see that the net effect of the BH at later times is to somewhat 
enhance the survival probability of binaries with binding energies 
of a few kTs (Eb < 16 kT), especially for binaries within the 
half-mass radius (see Fig. [SJ. This feature is especially apparent 
in a two-dimensional plot of the distribution of binaries as a func- 
tion of binding energy and radial distance from the cluster center, 
depicted in the bottom panels of Fig. [8] The rather strong corre- 
lati on between radius and b inding energy that is observed (e.g., 
see lGiersz & Spurzemll200ol Paperl) in systems without a central 
compact object is thus almost absent. The ratio of moderately hard 
(E b /kT < 16) to hard (64 < E b /kT < 128) binaries is greater by 
about a factor 2 at t ~ 20 t r h(0) in our runs with the presence of 
an IMBH (see top panel of Fig.|8j. The possible use of this ratio as 
a diagnostic to evidence the presence of an IMBH can be impaired 
if there is an initial segregation of primordial binaries depending 
on their semi-axis. Specifically if softer binaries are preferentially 
born in the core of the cluster then some of them could still survive 
after a few relaxation times, mimicking the signal due to an IMBH. 
However such a cluster would be characterized by a significantly 
smaller core radius with respect to a cluster truly hosting an IMBH. 

This picture, presented for a TV = 8192 Plummer model, is 
representative of all our set of simulations. We have shown in Sec|3| 
that the global evolution of the system is not sensitive on the details 
of the initial conditions. Similarly the evolution of the binary abun- 
dance at fixed relaxation times does not significantly depend on the 
number of particles TV (see Fig. [3} or on the details of the initial 
conditions: in Fig.|5|we depict the evolution of Nbi n /Nbi„(0) for a 
number of different simulations and the pattern is consistently sim- 
ilar. The same trend shown in Fig.|8|is also obtained starting from 
different initial conditions, as we have verified in the runs starting 
from Wo = 5 and Wo = 7 King models. All these results show 
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Figure 6. Location in the cluster of the disruption of hard binaries (mea- 
sured by KS regularization terminations) during the first 5 t rh (0) of a run 
with N = 8192 and 10% primordial binaries starting from a Plummer pro- 
file. In the left panel the radius is given in units of the local tidal radius for 
each disrupted binary (r< ), while in the right panel the radius is in units of 
the influence radius of the IMBH (r^). 
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Figure 7. Decomposition of the left panel of Fig.|6|into binaries disrupted 
inside (left) and outside (right) the influence radius n of the BH. The ori- 
gin of the bimodal distribution in Fig.|6]can be undestood in terms of an 
enhanced probability for a binary of experiencing a close encounter when 
located inside the influence sphere of the BH. 



that the evolutionary differences in runs with and without a central 
black hole are physical and not just due to random fluctuations. 



5 ESCAPERS 

The N-body systems in our simulations, despite being isolated, lose 
mass, i.e. some particles acquire a positive energy and are conse- 
quently removed from the system (if they are unbound and reach a 
distance of 20 th from the center of system). There are two main 
contributions to this process. One is diffusion of particles in the en- 
ergy space, that leads to the evaporation of particles with binding 
energy very close to (that translates into a small residual veloc- 
ity at infinity, which we call "ejection velocity"). The other is due 
to strong interactions of binaries with other stars and/or with the 
IMBH with the con sequent e jection of a high velocity particle (the 
slingshot effect, see lAarsethl2005l) . 

In our Plummer run with TV = 8192, / = 10% with tubh = 
0.014Af we find that a fraction of 3 ■ 10~ 3 of the stars of the cluster 
leave the system in a relaxation time and that w 1/3 of them have 
ejection velocities that exceed by at least five times the core veloc- 
ity dispersion. Among these 54 were binaries; in addition one triple 
was ejected from the system. A distribution of the ejection veloc- 
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Figure 8. Distribution of binding energies for binaries (upper panel) at 
time t = 24 t r ^ (0) for a simulation starting from a Plummer model with 
N = 8192, 10 % primordial binaries and m sff = 0.014M (solid line), 
compared to a similar simulation without the BH (dotted line). In the lower 
two panels we report the energy-radius distribution of binaries for the same 
simulations. The radii have been normalized to the Lagrangian radius con- 
taining 5 % of the mass in the simulation. The dashed line gives the position 
of the half-mass radius, at t = 24 i r ^(0). The relatively large survival at 
t = 24 t r h(0) of binaries within the half-mass radius with binding energy 
Ei, < 20kT is clearly visible in the bottom panel, and is also reflected in 
the higher values for the solid line at the left-hand side of the upper panel. 



ities, as measured in a series of 5 simulations with N = 8192, is 
shown in Fig. llQI 

These results, especially for the high velocity tail of the dis- 
tribution, need to be confirmed by more realistic simulations that 
include a mass spectrum and account for possible stellar collisions 
during close encounters, but they suggest that the interaction of 
IMBH with a significant population of binaries may efficiently lead 
to the production of a number of ultra fast stars, i.e stars with an 
ejection velocity of several hundreds km/s. 

The ejection velocity resulting from a binary-IMBH interac- 
tion is limited to be of the order of the orbital velocity of the 
binary members with respect to their center of mass (e.g., see 
lYu & Tremamell2003l) . This implies that for main sequence stars 
the maximum ejection velocity is of the order of the escape ve- 
locity from the star (of order of a few hundreds km/s) and this 
is approximately independent on the mass of the stars due to the 
(quasi) linear mass-radius relationship along the main sequence. 
Therefore only binaries with at least one highly compact member 
(a white dwarf, a neutron star or a stellar mass BH) can lead to the 
production of hypervelocity stars (i.e. stars with velocity exceeding 
lOOOfcm/s) in consequence of tidal shredding due to an IMBH. 
From the observational point of view this decreases the likelihood 
of a globular clusters origin for main sequence hypervelocity stars. 
The expected leading production channel is in this case an interac- 
tion with a supermass ive black hole to be produced, as initially pro- 
posed byjHm|]l988) and recen tly investigated by Yu & Tremaine 
j2003l) . lGinsburg & LoeblfcOOol) and lBaumgardt et alJ<2006l) . Un- 
fortunately the possibility of identifying observationally ejected 




10 15 

t/t rh (o) 

Figure 9. Like Fig.^but for a series of simulations (TV = 8192, / = 10%, 
rrXBH = 0.014) starting from King models with Wo = 5 (solid line; two 
different runs) and Wo = 7 (dotted line; two different runs). After an initial 
transient given by the different initial conditions, the binary destruction rate 
is very similar among all the runs and is quantitatively consistent with the 
one observed in Fig.|5]for an initial Plummer density profile. 



stars with velocities below 1000 km/s does not appear very promis- 
ing. In fact the several hundreds of km/s range of velocities is 
shared also by stars in the disk of our galaxy, so an accurate mea- 
sure of the proper motion would be required to trace back to the 
center of a globular cluster the trajectory of a candidate main se- 
quence ejected star. 



6 DISCUSSION AND CONCLUSIONS 

In this paper we have presented for the first time the results of direct 
N-body simulations of the evolution of a star cluster with a popu- 
lation of primordial binaries that harbors a central IMBH, with a 
mass of the order of 1 % of the total mass of the system. In order 
to begin to analyze the basic evolutionary processes we have re- 
stricted our simulations to isolated systems of stars of equal mass, 
neglecting at this stage stellar evolution and physical collisions. 

The environment around the BH turns out to be of great 
interest from a dynamical point of view. It represents a labora- 
tory where complex interactions between hierarchical systems take 
place. Around the BH we observe usually one or more stars tightly 
bound to it. With a significant population of primordial binaries, 
it frequently happens that a binary approaches the center of the 
system. Interactions between these subsystems are important not 
only from a dynamical point of view, but also because they form 
a factory to produce exotic objects, such as tight (X-ray emitting) 
binaries and high velocities escapers, which we observe in high 
numbers in our simulations. 

From a computational point of view the complex interactions 
between binaries and the BH have proved to be very challenging 
for current state of the art N-body algorithms. In our simulations 
we had to resort to several restarts of the computation to fine tune 
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Figure 10. Distribution of the escape velocity (in units of the core velocity 
dispersion at the escape time) as measured in a sample of 5 simulations 
with N = 8192, 10% primordial binaries, and m BH = 0.014. If applied 
to a typical star cluster (where <r c lOkm/s), this plot suggests that the 
combined action of an IMBH plus a population of primordial binaries can 
lead to the ejection of stars at several hundreds km/s. 



the integration parameters in order to achieve energy errors conser- 
vation below 5- 10 -3 at the end of the simulation and below 5- 10 -4 
per dynamical time. Clearly, new regularization methods have to be 
developed to provide both more accuracy and higher stability. This 
will allow to perform automatic runs, a necessary condition to em- 
ploy a more systematic study of systems with IMBH and binaries. 

We have shown that a BH with a mass of around 1 % of the 
cluster mass can produce most of the energy required to fuel the 
expansion of the star cluster, so that the evolution of the large scale 
structure of a star cluster with primordial binaries and a central BH 
is remarkably similar to that attained in the absence of binaries (see 
iBaumgardt et al]f2004al) . Even so, the evolution of the population 
of primordial binaries is strongly influenced by the BH: binaries 
are disrupted primarily within the influence sphere of the BH, but 
we note a significant lack of binaries tidally disrupted by the BH 
with respect to predictions based on the loss cone theory. These 
theoretical estimates only take into account interactions between 
the binaries and the BH starting from the density of binaries in the 
core. Instead we observe that most of the binaries are destroyed 
by three and four body encounters with other stars once they enter 
the influence sphere of the BH and before they have the chance to 
venture close enough to the BH to be tidally shredded. 

These three and four body encounters within the influence 
sphere of the BH often lead to the ejection of stars with proper 
motions easily exceeding 10 times the core velocity dispersion of 
the cluster. We observe events with velocities up to 50 times the 
core velocity dispersion, that is of the order of lOkm/s for a typ- 
ical globular clusters. Binaries within the influence sphere of an 
IMBH can therefore originate a significant population of high ve- 
locity stars ejected from the system, whose presence could be de- 
tected observationally. Tidal shredding of binaries with main se- 
quence components limits the ejection velocity to a few hundreds 



km/s, but higher velocities can be obtained if at least one mem- 
ber of the binary is a compact object (white dwarf, neutron star, or 
stellar mass BH). One related intriguing, although highly specula- 
tive possibility is indeed the dynamical origin of a w lOOOfcra/s 
pulsar candidate from the glo bular cluster remn ant IRS13, which is 
considered to host an IMBH <Wang et all2006l) . 

In principle, an indirect evidence for the presence of an IMBH 
in an old globular cluster could be obtained by observing the distri- 
bution of binaries: observational detection of an unexpectedly large 
fraction of only moderately hard binaries in the core could be used 
as circumstantial evidence for the presence of an IMBH in old star 
clusters, i.e. in those systems whose age would imply a depletion in 
the central region of binaries with binding energy below « 15 kT. 
With our runs we have shown that these binaries are able to sur- 
vive in significant numbers when an IMBH is present. For a typical 
globular cluster this would imply to be able to measure in the core 
the ratio of binaries with semi-major axis of about a > 0.5 AU to 
those with a ~ 0.05 AU, that appears to be particularly sensitive to 
the presence of an IMBH (see Fig.|8j. This relative measure has the 
advantage of being relatively independent from the precise value of 
the primordial binary ratio (clearly unknown for observed cluster), 
but it is extremely challenging for the present observational tech- 
niques of binary detection. This would in fact mean to be able not 
only to detect a significant number of binaries in the dense core en- 
vironment of a globular cluster, but also to accurately measure the 
semi-major axis of these binaries down to a fraction of AU. In ad- 
dition, as we have mentioned in Sec. 4, if softer primordial binaries 
are preferentially born in the center of the cluster, this diagnostic 
fails. 

A more promising observational evidence for the presence of 
an IMBH in old, relaxed, globular clusters can instead be given 
by a simpler measure, i.e. the ratio of the core to half mass ra- 
dius. This quantity is set by the efficiency of energy generation 
in the core of the cluster and is robust with respect to changes in 
the details of the initial conditions. Systems made of single stars 
only have a very small r c /rh ratio (r c /r^ « 0.01), as no effi- 
cient energy production mechanism is available. When primordial 
binaries are present, r c /rh is up to almost one order of magnitude 
larger and the ratio depends only on the initial fraction of bina- 
ries, saturating at / > 10% (see Paperland Paperll). Accounting 
for a small N-dependence of the ratio, a typical old globular clus- 
ter with primordial b inaries is e xpected to have r c /rh < 0.05 
(see also lFreseau & Rasid l2006). Most of the observed globular 
clusters have i nstead m uch larger values for r c /rh, located in the 
range [0.1; 1] iFregeau & Rasiol f2006). If we were to assume that 
old globular clusters host an IMBH, then we would naturally obtain 
r c /rh ~ 0.3 which is in good agreement with the observations. 

Clearly, a quantitative comparison with observations would re- 
quire us to proceed beyond the simple models presented here. To 
start with, we would need to include a realistic number of stars to 
avoid possible biases introduced by the scaling of the ratio of the 
BH mass to that of single stars. In addition, the introduction of a 
realistic initial mass function is likely to modify the concentration 
of binaries in the center of the system due to mass segregation, and 
stellar evolution will also influence the distribution of binary bind- 
ing energy. However, the main effects presented in this paper are 
likely to be present, at least in qualitative ways, in more detailed 
realistic simulations. 
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APPENDIX A: NUMERICAL ACCURACY 

The Aarseth's NBODY6 code represent one of the few state of the 
art integrators for the direct integration of collisional N-body sys- 
tems. Its performance and accuracy have been extensively tested 
(see lAarsetbl2003l) . The code employs a series of automatic inter- 
nal accuracy tests at the level of individual timesteps and forces, 
treatment of hierarchical systems and, of course, conservation of 
the total integral of motions. If one of these tests fails, the code 
resumes the computation from a previously saved snapshot of the 
system with improved integration parameters, to try to obtain an 
accuracy within the level required. If this is not possible, the code 
eventually stops the computation, returning an error message. 

Typically, for a system with N ~ 10 4 , the code can guarantee 
relative errors in the conservation of the total energy at the level 
of 10 -8 per dynamical time. The introduction of a population of 
primordial binaries increases significantly the dynamical range of 
the computation, as the shortest dynamical timescale to be resolved 
may be of the order of days, compared a dynamical time of a few 
million of years for a typical globular cluster. This is reflected with 
a relative energy error of the order of 10 -6 per dynamical time 
(see Fig. lAH . that still guarantees an excellent conservation of the 
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Figure Al. Relative total energy conservation per dynamical time for 
N = 8192 in a run with 10% primordial binaries (left panel) and in a run 
with 10% primordial binaries and an IMBH with mass 1.4% of the total 
mass of the system (right panel). The maximum allowed error per dynami- 
cal time is set at 5 • 10 -4 ; all our runs are stopped if the total energy error 
increases above 5 • 10 — 3 . 



energy at the end of the simulation (the errors tend to have ran- 
dom signs so that the total error increases approximately with the 
square root of the time). The code has a similar performance also 
for runs that include single stars plus a IMBH (see lBaumgardt et"al] 
l2004jbl) . 

The introduction of an IMBH in simulations with primordial 
binaries introduces additional challenges for the code, as in this 
case several multiple interacting systems have to be treated in the 
sphere of influence of the IMBH. Despite the improvements intro- 
duced to increase the performance of the code (see Sec|2j, signifi- 
cantly greater errors with respect to the standard runs for the code 
happened during our simulations (see Fig. lAH . leading to a num- 
ber of restart of the run performed manually in order to fine tune 
the integration options, depending on the kind of error message re- 
ceived. 

To ensure reliable results we set a maximum allowed energy 
error per dynamical time of 5 • 10 -4 and we definitely stop a sim- 
ulation that reaches a total energy error of 5 • 10" 3 .This threshold 
value is typically reached after about 25 initial half mass relaxation 
times (thus at a time longer than the Hubble time for a typical glob- 
ular cluster). 

This paper has been typeset from a TgX/ I5TgX file prepared by the 
author. 



